library(MetBrewer)
library(ggplot2)
library(dplyr)
library(tidyverse)
library(ggridges)
library(gghighlight) # highlights on ggplot2
library(here)
library(patchwork)
library(vdemdata)
library(magrittr)
library(ggrepel) # pointing tags and labels on ggplot
library(ggpubr) # writing stats on plots
library(scales) # rescaling variables package
library(wesanderson)
library(modelsummary)
library(openxlsx)

# Before running this script, run the script "psdi_computation.R"

vdem <- vdem%>%
  filter(year>1969 & year<2020)

parsys_vdem <- left_join(psdi.data.clean,vdem, by = c("country_name", "year"))

## Reverting the scale of government and opposition WHERE 0 IS LOW AND 1 IS HIGH
## to have same direction as the main PSDI variable
parsys_vdem$v2ps_democracy_government_inverted <- 1-parsys_vdem$v2ps_democracy_government
parsys_vdem$v2ps_democracy_opposition_adj_inverted <- 1-parsys_vdem$v2ps_democracy_opposition_adj

### Figure 2. CASES CONTENT VALIDITY ####
## PSDI EDI ###
psdi_edi <- ggplot(parsys_vdem, 
                   aes(x = v2ps_democracy_adj, y = v2x_polyarchy)) + 
  geom_point()+
  xlab("Party-System Democracy Index")+
  ylab("Electoral Democracy Index (V-Dem, 2023)") +
  labs(title = "V-Dem")+
  geom_hline(yintercept = 0.5, color = "black", linetype = "dashed") +
  scale_y_continuous(expand = c(0,0), breaks=seq(0, 1, by = .1)) +
  scale_x_continuous(expand = c(0,0), breaks=seq(-1, 1, by = 0.1)) +
  stat_cor(method = "pearson", label.y = .05, label.x = 0.15, r.digits = 3, p.accuracy = 0.001)+
  theme_light()+
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) 

psdi_edi <- psdi_edi + stat_smooth(method = "lm",
                                   formula = y ~ poly(x, 2))
# Polity5
p5v2018 <- read_excel("p5v2018.xls")
p5v2018$country_name <- p5v2018$country
parsys_p5 <- left_join(psdi.data.clean,p5v2018, by = c("country_name", "year"))

psdi_p5 <- ggplot(parsys_p5, 
                  aes(x = v2ps_democracy_adj, y = polity2)) + 
  geom_point()+
  xlab("Party-System Democracy Index")+
  ylab("Polity Score (Polity5, 2018)") +
  labs(title = "Polity5")+
  geom_hline(yintercept = 0, color = "black", linetype = "dashed") +
  scale_y_continuous(expand = c(0,0), breaks=seq(-10, 10, by = 1)) +
  scale_x_continuous(expand = c(0,0), breaks=seq(-1, 1, by = 0.1)) +
  stat_cor(label.y = -9.5, label.x = 0.15, r.digits = 3, p.accuracy = 0.001)+
  theme_light()+
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) 

psdi_p5 <- psdi_p5 + stat_smooth(method = "lm",
                                 formula = y ~ poly(x, 2))
psdi_p5
# FreedomHouse 
fh23 <- read_excel("All_data_FIW_2013-2023.xlsx")
fh23$country_name <- fh23$`Country/Territory`
fh23$year <- fh23$Edition
parsys_fh <- left_join(psdi.data.clean,fh23, by = c("country_name", "year"))

psdi_fh <- ggplot(parsys_fh, 
                  aes(x = v2ps_democracy_adj, y = Total)) + 
  geom_point()+
  xlab("Party-System Democracy Index")+
  ylab("Freedom in the World (FreedomHouse, 2023)") +
  geom_hline(yintercept = 50, color = "black", linetype = "dashed") +
  scale_y_continuous(expand = c(0,0), breaks=seq(0, 100, by = 10)) +
  scale_x_continuous(expand = c(0,0), breaks=seq(-1, 1, by = 0.1)) +
  labs(title = "Freedom House")+
  stat_cor(label.y = 5, label.x = 0.15, r.digits = 3, p.accuracy = 0.001)+
  theme_light()+
  theme(legend.position = "bottom",
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) 

psdi_fh <- psdi_fh + stat_smooth(method = "lm",
                                 formula = y ~ poly(x, 2))
psdi_fh

tiff('PSDI_democracy.tiff', units="in", width=14, height=7, res=600, compression = 'lzw')
psdi_edi + psdi_p5 + psdi_fh
dev.off()

### Figure 3. PSDI & DEMOCRACY LEVELS ####
cases <- parsys_vdem %>%
  filter(country_name == "Germany"| country_name =="Mexico"
         |country_name == "India" |country_name == "Russia")

ld.psdi <- ggplot(cases,
                  aes(year, 
                      group = factor(country_text_id.x), 
                      color = factor(country_text_id.x))) +
  geom_line(aes(y = v2ps_democracy_adj), color = "steelblue",  linetype="solid")+
  geom_line(aes(y = v2ps_democracy_codehigh_adj), color = "steelblue", linetype="twodash")+
  geom_line(aes(y = v2ps_democracy_codelow_adj), color = "steelblue", linetype="twodash")+
  geom_line(aes(y = v2x_polyarchy), color = "firebrick", linetype="solid")+
  gghighlight(cases$country_name =="Germany", 
              unhighlighted_colour = alpha("lightblue", 0),
              use_direct_label = F,
              label_key = country_name,
              label_params = list(size = 5))+ 
  xlab("Election Year")+ 
  ylab("Party-System Democracy Index") +
  scale_x_continuous(breaks=seq(1970, 2019, by = 10)) +
  labs(title = "Germany")+
  theme_light()+
  theme(legend.position = "bottom")

aut.psdi <-  ggplot(cases,
                    aes(year, 
                        group = factor(country_text_id.x), 
                        color = factor(country_text_id.x))
) +
  geom_line(aes(y = v2ps_democracy_adj), color = "steelblue",  linetype="solid")+
  geom_line(aes(y = v2ps_democracy_codehigh_adj), color = "steelblue", linetype="twodash")+
  geom_line(aes(y = v2ps_democracy_codelow_adj), color = "steelblue", linetype="twodash")+
  geom_line(aes(y = v2x_polyarchy), color = "firebrick", linetype="solid")+
  gghighlight(cases$country_name =="India", 
              unhighlighted_colour = alpha("lightblue", 0),
              use_direct_label = F,
              label_key = country_name,
              label_params = list(size = 5))+ 
  xlab("Election Year")+ 
  ylab("Party-System Democracy Index") +
  scale_x_continuous(breaks=seq(1970, 2019, by = 10)) +
  labs(title = "India")+
  theme_light()+
  theme(legend.position = "bottom")

demo.psdi <-  ggplot(cases,
                     aes(year, 
                         group = factor(country_text_id.x), 
                         color = factor(country_text_id.x))
) +
  geom_line(aes(y = v2ps_democracy_adj, color = "PSDI"),  linetype="solid")+
  geom_line(aes(y = v2ps_democracy_codelow_adj), color = "steelblue", linetype="twodash")+
  geom_line(aes(y = v2ps_democracy_codehigh_adj), color = "steelblue", linetype="twodash")+
  geom_line(aes(y = v2x_polyarchy, color = "EDI"), linetype="solid")+
  scale_color_manual(name = "Legend", values = c("PSDI" = "steelblue", "EDI" = "firebrick"))+
  gghighlight(cases$country_name =="Mexico", 
              unhighlighted_colour = alpha("lightblue", 0),
              use_direct_label = F,
              label_key = country_name,
              label_params = list(size = 5))+ 
  xlab("Election Year")+ 
  ylab("Party-System Democracy Index") +
  scale_x_continuous(breaks=seq(1970, 2019, by = 10)) +
  labs(title = "Mexico")+
  theme_light()+
  theme(legend.position = "bottom")

ea.psdi <-  ggplot(cases,
                   aes(year, 
                       group = factor(country_text_id.x), 
                       color = factor(country_text_id.x))
) +
  
  geom_line(aes(y = v2ps_democracy_adj), color = "steelblue",  linetype="solid")+
  geom_line(aes(y = v2ps_democracy_codelow_adj), color = "steelblue", linetype="twodash")+
  geom_line(aes(y = v2ps_democracy_codehigh_adj), color = "steelblue", linetype="twodash")+
  geom_line(aes(y = v2x_polyarchy), color = "firebrick", linetype="solid")+
  gghighlight(cases$country_name =="Russia", 
              unhighlighted_colour = alpha("lightblue", 0),
              use_direct_label = F,
              label_key = country_name,
              label_params = list(size = 5))+ 
  xlab("Election Year")+ 
  ylab("Party-System Democracy Index") +
  scale_x_continuous(breaks=seq(1970, 2019, by = 10)) +
  labs(title = "Russia and former Soviet Union")+
  theme_light()+
  theme(legend.position = "bottom")


tiff('PSDI_4cases.tiff', units="in", width=8, height=8, res=600, compression = 'lzw')
(ld.psdi+aut.psdi)/(demo.psdi+ea.psdi)
dev.off()

### Figure 4. PSDI CONVERGENT VALIDITY ####

### Mainwaring Perez-Linan
DDLA <- read_excel("DDLA.XLSX")
DDLA$country_name <- DDLA$country
parsys_ddla <-left_join(psdi.data.clean, DDLA, by = c("country_name", "year"))

parsys_npr <- ggplot(parsys_ddla,
                     aes(x = v2ps_democracy_adj, y = npr_all)) + 
  geom_point() +
  xlab("Party-System Democracy Index")+
  ylab("Normative regime preferences all actors, 1944-2010 (Mainwaring & Perez-Linan, 2013)") +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed") +
  scale_y_continuous(breaks=seq(-1, 1, by = .1)) +
  scale_x_continuous(breaks=seq(-1, 1, by = 0.25)) +
  labs(title = "A.") +
  stat_cor(label.x = 0, label.y = -.8, r.digits = 3, p.accuracy = 0.001)+
  theme_light()+
  theme(legend.position = "bottom") 

parsys_npr <-parsys_npr + stat_smooth(method = "lm",
                                      formula = y ~ poly(x, 2))

### power social group by Marquardt
ve_analysis <- readRDS("ve_analysis.rds")
ve_analysis$country_name <- ve_analysis$countryname
parsys_ve <-left_join(psdi.data.clean, ve_analysis, by = c("country_name", "year"))
# Marquardt's variable name is "z" (valid.a$z)
socgroupKM.plot1 <- ggplot(parsys_ve,
                           aes(x = v2ps_democracy_adj, y = z)) + 
  geom_point() +
  xlab("Party-System Democracy Index")+
  ylab("Identity-based Exclusion (Marquardt, 2022)") +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed") +
  scale_y_continuous(breaks=seq(-4, 4, by = .5)) +
  scale_x_continuous(breaks=seq(-1, 1, by = 0.25)) +
  labs(title = "B.") +
  stat_cor(label.x = 0, label.y = -2.5, r.digits = 3, p.accuracy = 0.001)+
  theme_light()+
  theme(legend.position = "bottom") 

socgroupKM.plot.fit1 <-socgroupKM.plot1 + stat_smooth(method = "lm",
                                                      formula = y ~ poly(x, 2))

tiff('PSDI_convergent.tiff', units="in", width=12, height=6.5, res=600, compression = 'lzw')
parsys_npr+ socgroupKM.plot.fit1
dev.off()

## Figure 5. Construct Validation ####

## left right 

## Applying same methodology used for PSDI but on economic left-right 
## from V-Party
lrps.raw <- vall %>%   dplyr::select(country_id, country_text_id, 
                                     country_name, year,v2pagovsup.dum,seatshare.1,
                                     v2pariglef) %>%
  dplyr::group_by(country_id, country_text_id, 
                  country_name, year,v2pagovsup.dum) %>% 
  dplyr::summarise(
    wm_lr = weighted.mean(v2pariglef,seatshare.1),
    .groups = 'drop') %>%
  as.data.frame()
## LRPS ###
lr.data <- lrps.raw%>%
  spread(v2pagovsup.dum,wm_lr) 
# rescale governmment and opposition to 0 to 1
lr.data$v2ps_economic_opposition <- lr.data$`0`
lr.data$v2ps_economic_opposition_adj <- ifelse(is.na(lr.data$v2ps_economic_opposition), 0,
                                               lr.data$v2ps_economic_opposition)
lr.data$v2ps_economic_government <- lr.data$`1`
lr.data$lrps <- lr.data$v2ps_economic_government+lr.data$v2ps_economic_opposition_adj
lr.data$v2ps_economic <- rescale(lr.data$lrps, to=c(0,1))
# drop Lithuania because IRT model on this case creates severe outlier
lr.data <- lr.data[lr.data$country_name != "Lithuania", ]

parsys <- left_join(psdi.data.clean, lr.data,  by = c("country_name", "year"))
parsys_demeco <- ggplot(parsys,
                        aes(x = v2ps_economic, y = v2ps_democracy_adj )) + 
  geom_point() +
  xlab("Party-System Left-Right Index")+
  ylab("Party-System Democracy Index") +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed") +
  scale_y_continuous(breaks=seq(0, 1, by = .1)) +
  scale_x_continuous(limits = c(.20,1),breaks=seq(0.10, 1, by = 0.1)) +
  labs(title = "A. ") +
  stat_cor(label.x = .2, label.y = .8,r.digits = 3, p.accuracy = 0.001)+
  theme_light()+
  theme(legend.position = "bottom") 

parsys_demeco <-parsys_demeco + stat_smooth(method = "lm",
                                            formula = y ~ poly(x, 2))
## Party System Institutionalization
psi_base <- read_csv("psi_base.csv")
parsys_psi <- left_join(psdi.data.clean,psi_base, by = c("country_name", "year"))

psdi_psi <- ggplot(subset(parsys_psi, v2ps_democracy_adj>0),
                   aes(x = v2ps_democracy_adj, y = psi)) + 
  geom_point() +
  xlab("Party-System Democracy Index")+
  ylab("Party System Institutionalization (Kim, 2023)") +
  scale_x_continuous(breaks=seq(-1, 1, by = 0.25)) +
  labs(title= "B.") +
  stat_cor(method = "pearson", label.x = 0, label.y = -2, r.digits = 3, p.accuracy = 0.001)+
  theme_light()+
  theme(legend.position = "bottom") 

psdi_psi <-psdi_psi + stat_smooth(method = "lm",
                                  formula = y ~ poly(x, 2))

# DPI 2020
DPI2020 <- read_csv("data/other data sources/DPI2020/DPI2020.csv")
DPI2020$country_name <- DPI2020$countryname
psdi.data.clean$country_text_id <- psdi.data.clean$country_name
DPI2020$system <- as.character(DPI2020$system)

newdpi <- DPI2020 %>% mutate(across(system,~as.factor(.))) %>%
  mutate(across(system,.fns = list(value = ~ as.numeric(.))))

## 1: -999; 2: Semi Presidential ; 3: Parliamentary; 4: Presidential
psdi_dpi <- left_join(psdi.data.clean, newdpi, by = c("country_name", "year")) 

dpi.plot <- ggplot(subset(psdi_dpi, system_value>1 & v2ps_democracy_adj>0), 
                   aes(x = factor(system_value), y = v2ps_democracy_adj,
                       fill = factor(system_value)))+
  geom_boxplot()+
  geom_signif(comparisons = list(c("2", "3"),
                                 c("4", "3"),
                                 c("2", "4")), 
              map_signif_level=TRUE)+
  geom_jitter(color="black", size=0.6, alpha=0.3) +
  xlab("Database on Political Institutions (Scartascini et al., 2021)")+
  ylab("Party-System Democracy Index") +
  scale_y_continuous(breaks=seq(0, 1, by = .25)) +
  labs(title = "C. ") +
  stat_cor(label.x = -0.7, label.y = -0.90 )+
  scale_fill_manual(values = wes_palette("Darjeeling1", n = 3),
                    labels = c("Semi-Presidential", "Parliamentary", "Presidential")) +
  labs(fill = "System: ") +
  theme_light()+
  theme(legend.position = "bottom")

tiff('PSDI_construct.tiff', units="in", width=12, height=12, res=600, compression = 'lzw')
(parsys_demeco+psdi_psi)/
  dpi.plot+
  plot_layout(heights = c(2, 1))
dev.off()


